Mapping gravitational lensing of the CMB using local likelihoods 
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We present a new estimation method for mapping the gravitational lensing potential from observed 
CMB intensity and polarization fields. Our method uses Bayesian techniques to estimate the average 
curvature of the potential over small local regions. These local curvatures are then used to construct 
an estimate of a low pass filter of the gravitational potential. By utilizing Bayesian/likelihood 
methods one can easily overcome problems with missing and/or non-uniform pixels and problems 
with partial sky observations (E and B mode mixing, for example). Moreover, our methods are local 
in nature which allow us to easily model spatially varying beams and are highly parallelizable. We 
note that our estimates do not rely on the typical Taylor approximation which is used to construct 
estimates of the gravitational potential by Fourier coupling. We present our methodology with a 
flat sky simulation under nearly ideal experimental conditions with a noise level of 1 /ii^-arcmin 
for the temperature field, /xiT-arcmin for the polarization fields, with an instrumental beam full 
width at half maximum (FWHM) of 0.25 arcmin. 



I. INTRODUCTION 

Over the past decade the cosmic microwave back- 
ground (CMB) has emerged as a fundamental probe of 
cosmology and astrophysics. In addition to the primary 
fluctuations of the early Universe, the CMB contains 
signatures of the gravitational bending of CMB photon 
trajectories due to matter, called gravitational lensing. 
Mapping this gravitational lensing is important for a 
number of reasons including, but not limited to, under- 
standing cosmic structure, constraining cosmolo^cal pa- 
rameters [m, [in and detecting gravity waves [igLTitI. |21| . 
In this paper we present a local Bayesian estimate that 
can accurately map the gravitational lens in high resolu- 
tion, low noise measurements of the CMB temperature 
and polarization fields. 

There is extensive literature on estimating the lens- 
ing of the CMB (classic references includ e j^, [l3l , [27| ) 
and some recent observational detections [3, [26|. The 
current estimators in the literature can be loosely char- 
acterized into two types. The first type was initiated in 
p7j (see also [El, [23|) and utilizes quadratic combinations 
of the CMB and its gradient to infer lens structure. The 
optimal quadratic combinations were then discovered by 
[isl, [20| and are generally referred to as 'the quadratic 
estimator'. This is arguably the most popular estimate 
of the gravitational potential and uses a first order Taylor 
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approximation to establish mode coupling in the Fourier 
domain which can be estimated to recover the gravita- 
tional potential (real space analogs to these estimators 
can be found in [H, Q). The second type is an approx- 
imate global maximum likelihood estimate and was de- 
veloped in [i,[9|. 

Our method, in contrast, locally approximates a 
quadratic form for the gravitational potential and es- 
timates the coefficients locally using Bayesian meth- 
ods. The locally estimated coefficients are then globally 
stitched together to construct an estimate of a low pass 
filter of the gravitational potential. The local analysis 
allows us to avoid using the typical first order Taylor ex- 
pansion for the quadratic estimator and avoids the likeli- 
hood approximations used in global estimates. Moreover, 
we are able to easily handle missing pixels, problems with 
partial sky observations (E and B mode mixing, for ex- 
ample), and spatially varying or asymmetric beams. The 
motivation for developing this estimate stems, in part, 
from current speculation that likelihood methods will al- 
low superior mapping of the lensing structure (compared 
to the quadratic estimator) under low noise levels, and 
that global likelihood methods can be prohibitively com- 
putational intensive — indeed intractable — without signif- 
icant approximation. 

We illustrate our mapping methodology on a high res- 
olution, low noise simulation of the CMB temperature 
and polarization field on a 17° x 17° patch of the flat 
sky. This simulation is used throughout the paper to 
demonstrate findings and techniques. To get an overview 
of the performance of our method. Fig. [T] shows the es- 
timated potential (left) from the simulated lensed CMB 
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FIG. 1. Left: Estimated gravitational potential on a 17^ x 17^ patch of the simulated flat sky. Right: Input gravitational 
potential used in the simulation. See Section |l] and Appendix [A] for the simulation details. 



temperature and polarization field (observational noise 
levels are set at 1 /ii^-arcmin for the temperature field, 
-arcmin for the polarization fields, with a beam 
FWHM of 0.25 arcmin). The input gravitational poten- 
tial is shown in the right diagram in Fig. [TJ The details 
of the simulation procedure can be found in Appendix [Al 
It is clear from Fig. [1] that the mapping accurately traces 
the true, unknown gravitational potential. To get an idea 
of the noise of this reconstruction for different realizations 
of the CMB + noise we present Fig. [2] which shows the 
different estimates of the projected matter power spec- 
trum using the estimated projected mass — with the local 
likelihood approach — for 10 different CMB + noise re- 
alizations (dashed lines) while keeping the gravitational 
potential in Fig. [1] fixed. The blue curve shows the esti- 
mated projected mass power spectrum if one had access 
to the true gravitational potential used in our simula- 
tions. Finally we plot the theoretical ensemble average 
projected mass power spectrum in red to get an idea of 
the magnitude of the errors in the mass reconstruction. 



II. LOCAL MAXIMUM A POSTERIORI 
ESTIMATES OF SHEAR AND CONVERGENCE 

The CMB radiation in the flat sky limit can be ex- 
pressed in term of the Stokes parameters T, Q, which 
measure total intensity T(ic), and linear polarization 
Q{x) and U{x) with respect to some coordinate frame 
X = (x, y) G M^. Instead of directly observing T, Q, U we 
observe a remapping of the CMB due to the gravitational 
effect of intervening matter. This lensed CMB can be 
written T{x + V0(cc)), Q{x + V(/)(cc)) and U{x + VMx)) 
where (j) denotes the gravitational potential (see Q, for 
example). 

To describe our estimate of the gravitational poten- 
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FIG. 2. Plot of the projected mass power spectrum (red) 
along with the estimated power spectrum using the true, but 
unknown, projected mass (blue). The dashed lines correspond 
to different estimates of the power spectrum using the esti- 
mated projected mass — with the local likelihood approach — 
for different CMB realizations but the same lensing potential 
realization. See Section |I] and Appendix lAl for the simulation 
details. 

tial, (/), first consider a small circular observation patch 
with diameter 5 in the flat sky centered at some point 
ico, denoted M5{xq) C . Over this small region we 
decompose (j) into an overall local quadratic fit and error 
term 

(/> = + e 

where is a local quadratic approximation of the po- 
tential (j) with error term e = (j) — . In what follows 
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we estimate ^ denoted ^ and associate this estimate 
with the neighborhood midpoint Xq. Then we repeat 
this procedure for other local midpoints Xq throughout 
the observation window. After a shrinkage adjustment 
is made to the local estimates we show, in Section IIID| 
how to stitch the estimates together to produce the final 
estimated potential (j) shown in Fig. [TJ 

Notice that as ^ ^ the expected magnitude of the 
error e approaches zero. This has the effect of improving 
the following Taylor approximation 



T{x + V0(x)) = T{x) + Ve{x) • VT(x) 



(1) 



for X G Msixo)^ where we use the notation x = 
X + Vq^ix) (with a similar Taylor expansion for both 
Q{x + V0(x)) and U{x + V0(x))). Notice that x de- 
pends not only on x but also the unknown coefficients 
of the quadratic term q^ . We brieffy mention that these 
are related to the convergence n and shear 7 = 7^ + ^72 
of the gravitational lens by 



K^-{qt + qt,)/2 

71 « -{€. 
72 « -qty 



yy' 

'/i)/2 



using the shear notation given in j27|. Now when 5 is 
sufficiently small we can truncate the expansion in ([T]) to 
get 



T{x + V(P{x)) 
q{x + V4>lx)) 
U{x + V(/)(a;)) 



T{x + Vq'^{x)) 
Q{x + Vq*{x)) 
U{x + Vq'^{x)) 



(2) 



on the local neighborhood A/5(iCo). By regarding q^ as 
unknown we can use the right hand side of (|2]) to de- 
velop a likelihood for estimating the coefficients of q'^ . 
Nominally q^ has 6 unknown coefficients for which to 
estimate. However, we can ignore the linear terms in 
q^ since the CMB temperature and the polarization are 
statistically invariant under the resulting translation in 
Vq^. Therefore, one can write q^ as ci(x — xo)^/2 + 
C2{x — xo){y — yo) + cs{y — ?/o)^/2 for unknown coeffi- 
cients ci = q^^, C2 = qty, C3 = q^y. 

An important probe of gravitational lensing from the 
CMB polarization is the creation of a curl-like B mode 
from the lensing [3, [HI- We remark that a local 
quadratic approximation in ([2]) still has the power to de- 
tect this B mode power so that the local procedure is not 
blind to this information source. To see this notice that 
a quadratic lensing potential remaps the coordinates by 



X 



X 



^xx ^xy 
^xy ^yy 



{x - Xo). 



If we assume the original polarization {Q{x),U{x)) is 
curl free then the lensed polarization has curl given by 

curl(Q(;r), U{x)) = -2j2U^{x) + 71 [Q,{x) - Uy{x)] . 



Therefore the shear parameter 7, and not the conver- 
gence is what creates local B-mode power. The dom- 
inant source of information for B-mode power is in the 
cross correlation between the lensed Stokes parameters 
Q{x) and U{x). This agrees with [l3| that the E-B cross 
estimator provides optimal signal to noise under nearly 
ideal experimental conditions. 

We finish this section with a remark on the accuracy 
of the Taylor approximation ([!]). As the the signal to 
noise ratio increases and the pixel resolution improves 
one can shrink the local neighborhood JVs{xo) so the 
term e becomes smaller (which improves the Taylor ap- 
proximation). However, as (5 ^ 0, the fields T, Q and U 
become nearly linear and one may expect some loss of 
information from the shrinking power in T, Q and U at 
frequencies with wavelengths smaller than the neighbor- 
hood JVs{xo). It therefore may be statistically advanta- 
geous to artificially increase the neighborhood size while 
simultaneously increasing the order of the local polyno- 
mial fit q^. Then, instead of recording the full polyno- 
mial fit at each midpoint Xq, one can retain the second 
order derivatives q^xi^o), qty{xo), q^{xo) for estimates of 

and 7. It is yet to be seen, however, what S and what 
polynomial order will be optimal for a given noise and 
resolution level. In Section IIIII we present an informa- 
tion metric for choosing the neighborhood size S for the 
simulation specifics and for a quadratic polynomial q^. 



A. The local posterior 

Using the Gaussian approximation of the CMB along 
with the quadratic potential approximation given by ([2]) 
we describe how to construct the likelihood as a func- 
tion of the unknown quadratic coefficients in q^. Let 
Xi, . . . ,Xn denote the observation locations of the CMB 
within the local neighborhood J^s{xo) centered at Xq. 
Using approximation ([2]), the CMB observables in this 
local neighborhood are well modeled by white noise cor- 
ruption of a convolved (by the beam) lensed intensity 
and polarization field. Let t, gr, u denote n- vectors of ob- 
served CMB values at the corresponding pixel locations 
in JVs{xo) for the intensity T and Stokes parameters Q, /7, 
respectively. Using Gaussianity of the full vector of CMB 
observables, z = {t^ ^q^ ,u^y , the log likelihood (up to a 
constant) as a function of the quadratic fit q^ can be 
written 

Ciq'^lz) = (E^^ ^N)~^z-^ Indet (E^, + TV) 

(3) 

where E^^ + TV is the covariance matrix of the ob- 
servation vector z (we use the subscript to emphasize 
the dependence on the unknown quadratic q^), N = 
diag(cr|./, cTg/, cr^/) is the noise covariance structure and 
/ is the nxn identity matrix. Notice that the noise struc- 
ture does not depend on the unknown quadratic q^. In 
the next section we will derive the exact form of the prior 
distribution on q^, denoted 7r{q'^), but briefiy mention 
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that the posterior distribution on ^ which we maximize 
to estimate ^ is 

oce^(^'l"V(g^). (4) 

The entries of E^^ + N contain the covariances 

(^fe^j)cMB' (^fe^i)cMB' (^fe^i)cMB ^ross covari- 

ances among t^q^u (we use tk to denote the k^^ entry 
of t, for example). Let cp denote the instrumental beam 
and (TT^o-Q^au denote the noise standard deviations of 
T^Q^U so that, for example, the k^^ entry of t is modeled 
as 



tk= (fxif{x)T{xk-x)^(7Tnk 



(5) 



where the n/e's are independent standard Gaussian ran- 
dom variables, Xk Vq^{xk) and x = x + Vq^{x). 
Note that this is an approximate model for tk based on 
(|2j). In actuality, the k^^ temperature measurement is 

]^2(P'X (f{x)T{Xk - X + V^(x/e -X))- 

earity of \/q^ allows us to write Xk - 
const ant +x/e — x on the small neighborhood A/i(xo). Un- 
der the assumption of zero B mode, the spectral densities 
associated with Q, U can be written 



f CTT^/e, but the lin- 
X + V(j){xk — x) 



^ = Cf cos(2(/p^) sm{2^i) 



(6) 
(7) 
(8) 



where tan((^^) = £2/^1 and i = (-^1,-^2) ^ I^^- Since one 
can write x + Vq^{x) = Mcc where the M is a 2 real ma- 
trix, the sheared Stokes parameters T(x), Q{x) and U {x) 
are stationary random fields with spectral densities given 
by C^_i^detM-\C^_i^detM-i and C^_i^detM-\ 
respectively. After adjusting for the beam (which is ap- 
plied after lensing) the covariance between the observa- 
tions in t can be written 



{tktj) 



3/CMB 



(27r)2 



'^^ ^' detM' 



(9) 

The computations are similar to complete the entries of 
covariance matrix + N. At face value the above in- 
tegral seems too computationally intensive for every pair 
Xk — Xj. Moreover, to apply Newton type algorithms 
for maximizing the posterior (|4]) one needs to compute 
the derivatives of (^fe^j)cMB with respect to elements of 
M. In Appendix [Bl we show that some of these compu- 
tational challenges can be overcome by utilizing a single 
FFT to quickly compute the above integral for sufficient 
resolution in the argument Xk — Xj to recover {tktj)Qy^-Q 
for all pairs k^j. 



B. Taylor truncation bias 

The quadratic function q^ is defined as the best least 
square fit of (j) over the neighborhood JVsixo). The resid- 
ual e = (j) — q^^ defined over A/'^(xo), is nonstationary 



and will therefore not have a spectral density that diag- 
onalizes the covariance structure. However, stationarity 
is a good approximation for order of magnitude calcu- 
lations on the truncation error in ([T]). We approximate 
the spectral density of e as an attenuated version of Cf 
by arguing that the quadratic fit effectively removes the 
spectral power at wavelengths greater than 26. Reason- 
ing similarly we expect the quadratic fit to have negligi- 
ble impact on the spectral power at wavelengths smaller 
than S. By assuming the spectral power grows linearly 
in the intermediary spectral range, from zero at ^ = tt/^ 
to C^^i^ at £ = 27r/(5, we obtain an approximate model 
for the spectral density of e 

C|«min{l,[^|£|-l]^}V 

where x'^ denotes the positive part of the real number 
X. Notice that the attenuation happens on the realiza- 
tions of (/), hence requiring the square on the low pass 
filter in the spectral density. This implies that the sec- 
ond term in the Taylor expansion ([T]) has approximate 
spectral density 



J 



detM 



In our simulation we use a neighborhood diameter 
of (5 = 0.006 radians (20.6 acrmin). This diameter 
was chosen using the information criterion developed 
in Section IIIII The corresponding approximate rms of 
VT{x) • Ve{x) is ~ 2.3 jiK with an order of magnitude 
reduction for the polarization field. Brute force simu- 

/ \ 
lation of ( Mean {T{xk + V0(x/e)) - T{xk)Y ) 

yields a value closer to ~ 3.6 /iiC, suggesting a reasonable 
stationary approximation to e. These approximations 
show that the polarization truncation error is smaller (by 
an order of magnitude) than the simulation noise level 
-arcmin. However, the temperature truncation 
error is greater than the temperature noise level 1 jj.K. 
A consequence is that the likelihood explains the addi- 
tional high frequency power in the observations (from 
the error term) by adjusting the estimate of q^ to artifi- 
cially magnify the convergence estimates. Indeed, this 
bias seems relatively constant and can be clearly seen in 
Fig.[3]in the top blue points. Each blue point corresponds 
to a local neighborhood: the x-coordinate representing 
the true V^g"^ associated with that neighborhood; the y- 
coordinate representing the estimated local value shifted 
up by 0.2, i.e. V^g^ + 0.2. The bias of nearly - 0.1 
above the top dashed blue line y = x + 0.2, shows the 
effect of the additional high frequency power of the er- 
ror term VT{x) • Ve{x). To adjust this, we subtract the 
overall mean of the local estimates, reasoning that the 
observation window is large enough at 17° x 17° so that 
the overall mean of the true values qi^^ , qi^y , is close to 
zero. For smaller observation windows it may be possible 
to estimate an overall quadratic fit to correct for this bias 
but we do not investigate that here. 



5 



Local estimated vrs truth 



0.5 



0.4 



0.3 



0.2 



0) 

I 0.1 



^ 



-0.1 



-0.2 



-0.3 



-0.4 



Raw estimates 




Aftei-- gradient- fit;. 



-. • ■■^^tT'^t'^^^ ■•• - After siirinkage correction 
./ V' ■ ■ • > . • •• 



-0.2 -0.1 0.1 0.2 

Simulation truth 



FIG. 3. Estimated values of V^(/)(cco), for each local neigh- 
borhood midpoint cco, plotted against the simulation truth at 
different stages of the algorithm. The blue points correspond 
to the raw estimates at each local neighborhood; The green 
points to the estimates after fitting a gravitational potential; 
The red points after a shrinkage correction. The y coordinates 
of the blue points are shifted up by 0.2 and the red points are 
shifted down by 0.2 to fit on the same diagram. See Sections 
III Bl and III Dl for discussion. 



C. The prior 7r((?^) 

The stationary approximation for e also yields an ap- 
proximation for the the prior distribution of the local 
quadratic fit q'^ using the identity = (j) — e. Since e 
is well modeled by a high pass filter of ^, the quadratic 
function q^ can be modeled by the corresponding low 
pass filter 



{£) 



J^s{xo)^ where = minjl, 

which has spectral density minjl, 



over X 

I + 



^l^lj I Cf. Therefore a natural candidate for the prior 

on the coefficients of q^ are the random variables %^^q^^. 
which are mean zero and Gaussian with variances ob- 
tained by the corresponding spectral moments of 
This prior is used on each local neighborhood Msixo) to 
derive the local maximum a posteriori estimate. For the 
simulation used in this paper, the neighborhood width 
was set to 5 = 0.006 radians (20.6 arcmin) which gives 
prior variances 0.0023, 0.0008, 0.0023 for g^^, and g^^, 
respectively (the only nonzero cross covariance is between 
(7^^ and q^y and is 0.0008). 



D. Reconstructing cj) from 

When observing the full sky, the estimates of ti will 
allow one to recover the gravitational potential (j) by 
solving the poisson equation V^(/) = —2k. (up to a con- 
stant). With partial sky observations, however, the shear 
is needed to break ambiguity corresponding to different 
boundary conditions. We do this in two stages, first using 
^xx^ Qxy Qyy (regarded as functions of the local neigh- 
borhood midpoint xq) to recover the estimated displace- 
ment field {(j)x^(l)y)^ then using this displacement field to 
recover the estimated potential (j). To handle this, we 

6oj as mini- 



adopt the method of p^, |2J and define 
mizers of functionals Fi and F2 defined as 



Fi{(l)x) = J dxdy [(0^^ - q^^)^ + ((/)3,^ - g^) 



F2{cl)y) = / dxdy\{cl),y - qtyf + {^yy - q^y) 



In particular, (j)x satisfies Fi{(j)x) = min^^ Fi(0a^) and 
similarly for ^^y- See ^|] for details of the minimiza- 
tion algorithm. Now we use the estimated displacements 
{(pxi^Py) to define the estimated potential (j) as the mini- 
mizer of the functional F3 defined as 



F3((/)) = / dxdy 



■ C^xf + {^y - C^yf 



The minimization is needed to account for the fact that 
our estimates are noisy versions of the truth and therefore 
may not correspond to an integral vector field for which a 
potential exists. A consequence is that the estimate (j) is 
'shrunk' towards zero when the algorithm fits a gradient 
to a vector field which may have non vanishing curl. This 
shrinking can be seen in Fig. [3] looking at the scatter plot 
of green points. These points show (V^0(xo), V^^(iCo)) 
for each local neighborhood A/'^(xo). One can clearly see 
the shrinkage effect by noticing the slope of the trend in 
the green points is less than one. We undo this shrink- 
age effect by multiplying ^ by a factor that undoes this 
shrinkage. The multiplication factor, denoted c, is de- 
termined by matching the variance of the raw estimates 
V'^q^ with cV'^cj). The result of this correction factor is 
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FIG. 4. The right diagram shows (j)xj where = {(^x, (j^y) is the true gravitational displacement field used in the simulation. 
The left diagram shows the estimate (j)x which is derived from the local quadratic estimates using the methodology described 
in Section Til Dl 



seen in the scatter plot of the red points, in Fig. [3l which 
show the local convergence estimates versus truth after 
the correction factor (V^0, cV^^ — 0.2). 

The estimated (after correcting for the shrinkage) 
along with the true displacement (j)x (used in the simula- 
tion) are shown in Fig. HI The estimated cj) along with the 
true gravitational potential (j) are shown in Fig. [H These 
two figures demonstrate accurate reconstruction of both 
the gravitational potential and the displacement field. In 
addition, by differentiating the estimated potential, ^, 
one obtains smoothed estimates of convergence and shear 
(smoothed from the fitting of (j)). In Fig. [5] we plot the 
estimate (j)xy (which corresponds to minus the imaginary 
part of the shear 7), along with (p^y (bottom right) and 
the low pass filter (j)^^y (bottom left) defined in Section 

III C[ Notice that the estimate (j)xy tracks the derivatives 
of the low pass filter 0^^^, whereas the additional high 

frequency in (p^y is not accurately estimated from (p^y 
This is presumably due to the local fitting of a quadratic 
potential over the neighborhoods Msixo). 



III. NEIGHBORHOOD SIZE AND STRUCTURE 

We define the following measure of information which 
is used as a metric for choosing the width of the neigh- 
borhood and other parameters of our estimation method: 

Information for = 

variance of the prior on 
expected variance of the posterior on q^ 

The above information metric is essentially a measure 
of signal to noise ratio (squared). The variance of the 



prior corresponds to the squared magnitude of the sig- 
nal, whereas the expected variance of the posterior is a 
proxy for the squared magnitude of the noise. We use 
simulations to estimate this information (while using the 
hessian of the posterior density at q^ to approximate pos- 
terior variance) and use it for guidance when choosing the 
tuning parameters for our estimation algorithm. Note: 
we avoided a lengthy and rigorous simulation study to 
choose global optimal tuning parameters, opting for a 
less rigorous simulation study which yields, potentially, 
sub-optimal but reasonable algorithmic parameters. 

The main parameter that needs tuning is the local 
neighborhood size S. Notice that our information mea- 
sure attempts to balance two competing quantities when 
choosing a neighborhood size, the larger the neighbor- 
hood the smaller the signal q^ (from the low pass filter). 
On the other hand, larger neighborhoods correspond to 
more data when the resolution is fixed. Using this met- 
ric, 6 = 0.006 radians (20.6 arcmin) emerges as a good 
neighborhood size when the beam FWHM is 0.25 arcmin 
and the noise levels are \/2 and 1 /ii^-arcmin pixels for 
Q, U and T, respectively. 

Due to computational limitations associated with 
larger neighborhoods we found it necessary to down- 
sample the local neighborhoods by discarding pixels. Us- 
ing the information metric we were able to isolate that 
randomly sampling the pixels seemed preferable to evenly 
downsampling to a courser grid. Moreover, we found that 
using different randomly selected pixels for T, Q and U 
was preferable to using the same random pixels for all 
the Stokes fields. Therefore, for each local neighborhood 
we selected 300 random pixels in JVs{xo) for the T ob- 
servations, then randomly selected 300 pixels from those 
remaining for Q and finally 300 pixels from the remain- 
ing unselected pixels for U (allowing overlaps when the 
local neighborhood size had fewer than 900 pixels). 
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FIG. 5. The top diagram shows the estimate of (\)xy (which corresponds to minus the imaginary part of the shear 7) where 
denotes the gravitational potential. The bottom two diagrams show the simulation truth: bottom left shows (^^y where 

(^^ denotes the low pass filter = min|l, — Q-QQ^ (see Section [ll CI for a discussion); bottom right shows 

(l)xy Notice that the estimate of (l)xy tracks the low pass filter 0^^^ and does not have the high frequency behavior seen in the 
simulation truth Sxv. 



IV. DISCUSSION 



We have demonstrated the feasibility of using a local 
Bayesian estimate to accurately map the gravitational 
potential and displacement fields under low noise, small 
beam experimental conditions. The motivation for de- 
veloping this estimate stems, in part, from speculation 
that likelihood methods will allow superior mapping of 
the lensing structure (compared to the quadratic estima- 
tor) under low noise levels. The main difference between 
the global estimates of [H, [9[ and the local estimate pre- 
sented here is the nature of the likelihood approximation. 
In jl, ^ the global likelihood is defined as a functional 
on the unknown gravitational potential and approxima- 
tions are made to this functional. Our method, in con- 
trast, uses a nearly exact likelihood — exact up to approx- 



imation (jBip in Appendix [B] — but under a local model- 
ing approximation that assumes a quadratic (j). One ad- 
vantage is the added precision available to model instru- 
mental and foreground characteristics. For example, the 
local analysis models the beam convolved CMB rather 
than the deconvolved CMB. Deconvolution induces spa- 
tial correlation in the additive instrumental noise which 
is potentially nonstationary if the beam spatially varies. 
Since this noise is not invariant under warping it com- 
plicates the global likelihood. Another advantage is that 
the local estimates are relatively easy to implement and 
parallelize. In addition, the local estimate automatically 
uses the highest signal to noise combinations of Q, /7 and 
T so there is no need to re-derive the optimal quadratic 
combinations for different experimental conditions. 

The local analysis is not free from disadvantages how- 
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ever. A global analysis is presumably much better suited 
for estimating long wavelengths in the gravitational po- 
tential and wavelengths that are shorter than the lo- 
cal neighborhood size. Moreover, since our estimates 
are defined implicitly — as the maximum of the posterior 
density — it is difficult to derive expected error magni- 
tudes. However, the results presented here show that 
under some experimental conditions the advantages over- 
come the disadvantages. Moreover our local estimate 
uses an approximation that is inherently different from 
the Taylor approximation used to derive the quadratic 
estimator. This leaves open the possibility that the local 
estimate may have different bias and error characteristics 
which could compliment the quadratic estimator, rather 
than replace it. 



{wlwir) = (wiwf) = (wlwf) = 0, 

W^i = and Wl^ = Wf. 

In our simulation, the above sums — equations 
(|Aip . (|A2[) and (|A3[) — are taken over frequencies i G 
{^k : k G {-7V/2,...,iV/2 - l}^} where L = 0.2967 
radians so that T will be periodic on [— L/2, 1//2]^. The 
limit N = L/ is chosen to match the resolution in pixel 
space, denoted A^^, so that FFT can be used to compute 
the sums (|A1[) . (|A2[) and (|A3[) which, after simplification, 
becomes 



t'= L 
2n 



(A5) 



QUA,) « Yl Zl^—e^^^''-^/'' cos(2^,)^/c|^ (A6) 



Appendix A: Simulation details 



— '^ L 



The fiducial cosmology used for the simulations is 
based on a fiat, power law ACDM cosmological model, 
with baryon density = 0.044; cold dark matter 
density l^cdm = 0.21; cosmological constant density 
Qa = 0.74; Hubble parameter h = 0.71 in units of 
100 kms~^ Mpc~^; primordial scalar fluctuation ampli- 
tude As{k = 0.002 Mpc"^) = 2.45 x 10"^; scalar spectral 
index ns{k = 0.002 Mpc~^) = 0.96; primordial helium 
abundance Yp = 0.24; and reionization optical depth 
Tr = 0.088. The CAMB code is used to generate the 
theoretical power spectra [18]. 

We start by simulating maps of the unlensed CMB 
Stokes parameters T^Q^U . The following Riemann sum 
approximation is used for the random fields T^Q^U 



T{x) « 



2 Ax t 



27r 




where ipi = tan~-^(-£2/-^i); A£i,A-^2 are the frequency 
spacing in the two coordinate directions; for each £, 
Zj and Zf are mean zero complex Gaussian random 
variables such that {ZjZj*) = {ZfZf*) = 5f_ 



c 



Zj" and Zf^ 



{Z'^.Zf) 



To enforce the proper cross correlation between Zj and 
zf we set 



r 7^ 1 


1 






"71 





(A4) 



where p = 



and for each W, 



wl 



are 



for each j G {7V/2, . . . , Ar/2 
over k G {-7V/2, . . . , Ar/2 ■ 
ues {T^3^.)\.^^_^i^_^i^_^y 
be simulated by a two dimensional FFT of the matrix 



1}^ where the sums range 
1}^. The matrix of val- 
for example, can then 



27r i27rfc-j7Ar 



J ke{- 



. The iden- 

Ar/2,...,Ar/2-i}2 

titles W\j^ = wl* and = are enforced using a 
two dimensional FFT of two N x N matrices with inde- 
pendent standard Gaussian random entries. 

Remark: Typically the above method suffers from 
an aliasing error when truncating to a finite sum in 
()Aip , ()A2p and ()A3p . We avoid any such complication by 
setting the power spectrum in Cj, Cf and to zero for 
all frequencies beyond \i\ = 6000. We justify this trunca- 
tion since both diffusion damping and the beam FWHM 
of 0.25' combine to produce negligible amplitude in the 
CMB Stokes parameters at frequencies \£\ > 6000 com- 
pared to the noise level. 

Remark: Since the full sky Stokes parameters T^Q^U 
are defined on the sphere, the theoretical power spectrum 
C?, CY are only defined on integers I. Our 



for Cj 



mean zero complex Gaussian random variables such that 



flat sky approximation is obtained by extending Cj^ 
and to £ G by rounding the magnitude |£| to the 
nearest integer. See Appendix C in [11] for a derivation 
of this flat sky approximation. 

To get a realization of the lensed CMB Stokes parame- 
ters T, Q, /7 we use the above method to generate a high 
resolution simulation of T, Q, and the gravitational po- 
tential (/) on a 17'' X 17'' patch of the flat sky with 0.25 
arcmin pixels. The lensing operation is performed by 
taking the numerical gradient of 0, then using linear in- 
terpolation to obtain the values T{x + V(j){x))^Q{x + 
V^(cc)), /7(ic + V(j){x)). The beam effect is then per- 
formed in Fourier space using FFT of the lensed fields. 
Finally, we down-sample the lensed fields, every 4*^ pixel, 
to obtain the desired arcmin pixel resolution for the sim- 
ulation output. 



9 



Appendix B: Newton's method for maximizing the 
local posterior 

In this section we discuss our numerical procedure for 
maximizing the local posterior given by (j4]). We remark 
that calculations need to be fast since they will be per- 
formed on each local neighborhood for which a shear and 
convergence estimate is required. We discuss how the 
FFT can be used to to compute the covariance matrix, 
denoted + in Section III A( and the corresponding 
derivatives with respect to the unknown coefficients of 

. We let M be the symmetric 2x2 matrix defined as 

^xy ' ^yy y 
matrix M is regarded as the unknown which will be esti- 
mated from the data in the local neighborhood JVs{xo). 
Let Ti)^M denote a sheared temperature field, convolved 
with a Gaussian beam (with standard deviation a^) so 
that 



M 



so that + V(7^(ic) = Mx. The 



(x) 



R2 



(fyT{Mx - My)e 



-|yr/(2a,^) 



(ag2^)- 



To compute the covariance matrix of the T observations 
t = (ti, . . . in M3{x{)) (see equation (j5j)) one needs 
to evaluate the following covariance function for a given 
test shear matrix M at all vector lags h = Xj — Xk 



-j 



{Tb,M{x - 



h)Tb,M{x))cMB 



detM 



Once the covariance matrix E^^ + is constructed us- 
ing the approximation ()Bip (and the analogous approxi- 
mations for Q^U and all cross correlations) the posterior 
is easily computed as p{q^\z) ex e^*^^ ^^^7r{q^) where C 
denotes the log likelihood ([3]) and tt is the prior distri- 
bution derived in Section III CI In principle, one can now 
simply use pre-existing minimization algorithms for max- 
imizing the posterior p{q^\z) with respect to q^. If one 
desires a more sophisticated Newton type algorithm for 
maximizing the posterior one often needs to compute the 
gradient and hessian of the posterior. Using the tech- 
niques of automatic differentiation (see [19], for exam- 
ple) one can easily compute such derivatives if one can 
compute the rates of change of the covariance E^^ with 
respect the the elements of M. 

We finish this Appendix by noticing that the FFT can 
be used to approximate the derivatives of E^^ with re- 
spect to the elements of the matrix M, denoted M^j for 
/c, j G {1,2}. For illustration we focus on the covariance 
structure of the temperature field T and mention that 
the extension to Q^U is similar. First notice that by 
transforming variables t = M~^i one gets 



dMkj 



d 



dH' 



dMkj J (27r)2 
dH' 



I 



^i{M£.').h-cjl\M£.'\^^T 
i{Mi')-h-al\M£'\'^ f^T 

i(M(!)-h-al\M(!^ 



d 



dMk 



All these calculations can be approximated using a FFT 
by noticing 



Cr,,„(jA,) 



E 

k 



A2 



^i27rk-j/N^- 



detM 
(Bl) 



where the sum ranges over k G {—A/2, . . . , A/2 — 1}^, 
Aa, is the pixel spacing, j e {-A/2, . . . , A/2 - 1}^, 
= 27r/L and L = NAx. Then to compute the co- 
variance between tj and tk we simply select the entry 
of the matrix [Cti^mU^x 
such that jAx 



Ar/2,...,Ar/2-l}2 
Xk (which was obtained by a sin- 



gle FFT). A similar technique can be used to compute all 
other covariance and cross-covariances among T, Q and 
U to construct the covariance matrix E^0. We remark 
that to speed up the computations we choose a smaller 
A then the one used in the simulations (A = 4096 in 
the simulation but A = 256 for the approximation of 
Cr.^ijAx)). 



Now d 



i{Mi')-h-(jl\Mi'\^ 



dMkj can be written as a 



sum ^^Ck(h)gk{M^l') so that by re-transforming vari- 
ables to i = Mi' one gets 



dM, 



(27r)2 



' ^detM 



The point is that the above integrals can now be approx- 
imated using FFT to approximate the matrix of values 

^^^^J^/r^^^"^^ • The same method ap- 

(^^^^k,3 J je{-A/'/2,...,Ar/2-l}2 

plies to approximate all higher order derivatives of the 
covariance matrix. These derivatives can then be used 
in a Newton type algorithm for finding the maximum 
a-posteriori estimates q^. 
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